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SUMMARY 

Finite element substructuring strategies to obtain numerical solutions 
for three typical problems of interest to the composites community, are pre- 
sented in this paper. The key issue common to these problems is the presence 
of singular or near singular stress fields. The regions prone to see steep 
stress gradients, are substructured with progressively refined meshes to study 
the local response simultaneously with the global response. The results from 
the select examples indicate that finite element substructuring methods are 
computationally effective for composite singularity mechanics. 


INTRODUCTION 

Finite element substructuring analysis has long been recognized as a via- 
ble alternative to analyze large scale structures. Efficient use of computer 
storage, reduction in computation time, and much less turnaround time are some 
of the advantages of finite element substructuring. Over the past few years, 
at the NASA Lewis Research Center, some of the solution capabilities available 
in MSC/NASTRAN have been utilized to analyze several problems related to com- 
posite structures (refs. 1 to 4). This research has focused on singular or 
near singular stress field behavior in composites. The approach has been to 
substructure the local regions with a progressively refined mesh and study the 
local response simultaneously with global response. The regions chosen for 
progressive refinement are near load application points, crack tip zones, and 
free edges. It is a well known fact that near-singular or singular stress 
fields exist in these regions. Experiences with the use of finite element sub- 
structuring solution strategies are described via three typical problems which 
represent current issues of interest in the composites community. However, 
the underlying solution principles are not limited to such problems only. 

The selected problems are (1) impact and toughness characterization of 
composites using Charpy's impact test specimen, (2) free-edge stress analysis 
of composite laminates, and (3) fracture toughness predictions of composites 
for individual and combined fracture of modes I, II, and III. The objective of 
this paper is to report on the effectiveness of the finite element substructur- 
ing concept in conjunction with the analysis of structural composites with 
embedded singularities. Another important objective of the present study is to 
see how well one could relate the global response characteristics to the con- 
stituent level properties (fiber/matrix micromechanics) using two- and 
three-dimensional finite element analysis. 



SELECT EXAMPLES 


Impact and Toughness Characterization of Composites 
Using Charpy's Impact Test Specimen 

Specimen geometry and finite element idealization . - The geometry of the 
standard Charpy test specimen (ASTM STD E23-7) along with the finite element 
idealization for smooth and notched specimens are shown in figure 1. Two 
finite element models of Charpy test specimens— one without a notch (smooth 
specimen) and one with a notch were selected for the present study. The 
material properties are assumed to be uniform, orthotropic, and obey a linear 
stress-strain law throughout the analysis. In addition, the specimen is 
assumed to be in a state of plane stress. The plane stress assumption is jus- 
tified from the physics of the problem. The width restraints at the notch-tip 
are negligible because of the very low value of the respective Poisson's ratio. 
For the present analysis, the plane stress assumption implies that the stresses 
are permitted to vary along the specimen length and through the thickness but 
not across the width. This reduces the stress calculations to three; two nor- 
mal and one shear. 

With these assumptions, plane stress finite elements can be used to model 
the Charpy test specimen. For the smooth specimen all the elements are quadri- 
lateral. The notched specimen is modeled with both triangular elements and 
quadrilateral elements. The triangular elements are used as transition ele- 
ments in the areas around the supports, the load application point, and the 
notch. These are the regions where maximum stress concentrations are expected 
to occur and therefore are modeled with a finer mesh. The boundary conditions 
prescribed are such that the node at the left support is constrained from x, 
y, and z displacements, and the node at the right support is constrained 
from y and z displacements. In addition, for the notched specimen, three 
nodes closest to the right support are constrained from displacement in the 
y direction. The specimens are subjected to an impulse loading. The form of 
impact is a triangular function with a peak load of 2000 lb occurring at 500 ys 
of a total contact time of 1000 ys. 

The details of the finite element models for the smooth and notched speci- 
mens are described in reference 1. 

Finite element analysis method . - The MSC/NASTRAN general purpose struc- 
tural analysis finite element computer program is used for the finite element 
analysis. The specific elements used are identified as CTRIA3 and CQUAD4. 

They are isoparametric elements. NASTRAN obtains the solution using a dis- 
placement formulation via rigid format solution sequence number 27. This solu- 
tion sequence employs a direct time integration scheme to obtain the transient 
response of a structure subjected to impact. 

The solution sequence number 27 of MSC/NASTRAN uses the integration algo- 
rithm based upon the Newmark Beta method (ref. 5). It provides stable results 
for the widest possible spectrum of practical problems without sacrificing 
either accuracy or efficiency. For complete details of the MSC/NASTRAN analy- 
sis, reference 6 should be consulted. 

Composite system analyzed . - Flexural (Charpy-type) test specimens made 
from three typical composite systems are analyzed. They are T-300/Epoxy, 
Kevlar/Epoxy, and S-Glass/Epoxy composites. The specimens are all 
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unidirectional composites with the fibers parallel to the length (x-axis, 
fig. 1 ) of the specimen. 

The plane stress-strain relationship (stiffness) coefficients required to 
input to NASTRAN are obtained by using the resident data-bank in the composite 
mechanics computer code ICAN (ref. 7). 


Two types of load conditions are used in obtaining the results. The first 
is a static loading where a force of 2000 lb is applied at the center of the 
specimen on the top surface. The second is an impulse loading. This modeled 
as a triangular pulse with a peak value of 2000 lb in the middle. The pulse 
is modeled to last for 1000 ps. The transient response is, however, obtained 
for three contact time periods i.e., 3 ps. 


Displacement and stress wave propagation . - The bulk wave and shear wave 
velocities are normally much higher compared to the flexural wave velocities. 

In order to capture the characteristics of propagation of these waves, a much 
smaller time step of integration (0.1 ps) is chosen. Figure 2 shows the 
dynamic displacement propagation in S-Glass/Epoxy, Kevlar/Epoxy, and T-300/ 
Epoxy specimens. Two bulk wave velocity parameters and one shear wave velocity 
parameter are defined below to aid the following discussion of the displacement 
wave propagation results: 



(1) 


( 2 ) 


(3) 


where p is the mass density of the material. 

The velocities are expressed in in./ps. The trend indicated by Cg22 f° r 
the three materials (the transverse shock wave travels fastest in S-Glass/Epoxy 
and slowest in Kevlar/Epoxy) is seen clearly in figure 2. A rough estimate of 
the normal wave velocity can be obtained by counting the number of elements 
that appear to be affected by the impact from figure 2. The normal wave veloc- 


ty estimates from the 
lacement propagation 

1 ps and 3 ps frames are shown 
after 1 ps is not shown in fig. 

below: 

2) 

(note that dis- 

Composite system 

Number 
1 ps 

of elements 
3 ps 

Velocity 
1 sec 3 sec 

Average 

S-Glass/Epoxy 

10 

25 

0.1313 

0.1094 

0.1204 

T-300/Epoxy 

7 

16 

.0788 

.0700 

.0744 

Kevlar/Epoxy 

5 

12 

.0657 

.0520 

.0589 
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These values are in close agreement with the values predicted using equa- 
tion (2) for Cq22 • The same trend is also seen for the 3 and 5 ps frames. 

It can be concluded that the initial shock travels with the bulk wave veloci- 
ty Cb 22 along the direction of impact. 


Once the normal shock reaches the bottom of the beam, the wavefront 
appears to be moving in the longitudinal direction forming a flexural wave. 
The velocities of the waves traveling in the longitudinal direction can be 
estimated with the same technique mentioned earlier. The following are the 
details for the frame after 13 ps (not shown in fig. 2). 


Composite system 


Number of elements 


Velocity, 

in./ps 


S-Glass/Epoxy 

T-300/Epoxy 

Kevlar/Epoxy 


11 

9 

7 


0.0510 

.0416 

.0320 


The above velocities have the same trend as the bulk shear wave velocity 
parameter C$i 2 • However, the waves propagating in the longitudinal direction 
appear to move significantly slower than the calculated shear wave velocity. 
This is probably due to the coupling between the flexural wave and the shear 
wave. The flexural wave velocity is significantly slower than the shear wave 
velocity. For example, the smooth S-Glass/Epoxy specimen has a flexural wave 
velocity Cf given by 


Cp = 2Pf = 0.0137 in./ps 


where P is the length between the supports and f is the first fundamental 
frequency. It is assumed that the wave number is 1 and the beam deflects into 
a half wave. 

Static and dynamic stress contours . - The static and dynamic stress con- 
tours for notched S-Glass/Epoxy specimen are shown in figure 3. They appear 
to be identical. The peak values for the dynamic load case, however, are 
slightly lower than those for the static case. Steep stress gradients are 
observed near the load application point and at the notch tip. It is seen that 
all the three stresses attain very high peaks, indicating severe local stress 
intensities near the notch-tip. As the stress allowables for <722 and 012 
are generally an order of magnitude lower than that for an, one can except a 
matrix initiated failure followed by fiber fractures at this location. 


Free-Edge Stress Analysis of Composite Laminates 

Specimen geometry and finite element idealization . - The geometry of the 
laminated plate under study and the finite element mesh employed for the pri- 
mary and superelement structures are shown in figure 4. 

One element is used through the thickness of each layer. There are 7 ele- 
ments across the plate width and 28 elements along the plate length. A total 
of 1365 brick elements are used to model the primary structure. The superele- 
ment contains 224 elements with progressively decreasing mesh size toward the 
free edge. With this type of arrangement, the center of the element nearest 
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to the edge is 1/512 in. from the edge. The boundary conditions imposed are 
such that along the edge x = 0 every node in the midplane is constrained from 
displacements in x, y-, and z-di rections . Also, the rest of the nodes along 
the edge x = 0 are constrained from displacement in the x-direction. Along 
the edge x = 28 in. a uniform stress of magnitude unity is prescribed. 

Cases investigated . - The types of laminates investigated in the present 
study are <±10) s through (±80) s at 10° intervals. The (±45) s laminate is also 
studied as a special case since it has been given substantial attention in the 
literature. The interply layers are assumed to be isotropic with E » 0.5 mpsi 
and v = 0.35. These are representative properties of an epoxy matrix. The 
thickness of the interply layer is taken as 0.000067 in. The plies are assumed 
to be homogeneous and orthotropic. The material properties chosen for the 
plies are typical of an AS-graphite fiber/epoxy composite system and are gener- 
ated using, ICAN (ref. 7). The details are described in reference 2. 

Free-edqe stress . - In order to study the three-dimensional stress behav- 
ior near the straight free edge, a small region close to edge is isolated and 
modeled as a substructure (superelement). This region is shown in figure 4. 
Typical results for a (±10) s AS-graphi te/epoxy laminate are shown in figure 5. 
The six stresses, which were obtained near the free-edge region, are summarized 
in figure 5. These results are plotted against a length parameter Xf along 
a central line perpendicular to the free edge. The length parameter Xf is 
related to the distance from the free edge Lf (fig. 4) by 


Xf = 1 - Lf (4) 

In the figure the dashed lines represent the stress computed at the center of 
(±9) ply elements. The solid lines represent the stress computed at the cen- 
ter of the interply layer between (+9) and (-9) layers. The stresses are nor- 
malized with respect to the applied stress o XX0 0 . The classical laminate 
analysis predictions for these normalized stresses are 

- l; ^- = 0 (5) 

o o 

XX® XX® 

Inplane stresses aw, a V y and ay y. - The inplane stresses a xx , o X y, 
and oyy approach classical laminate theory predictions in the interior or the 
ply. Close to the free edge, however, significant changes in these stresses 
are clearly seen. The transverse normal stress oyy in the +9 ply appears 
to be negligible. However, it does not vanish because of the presence of the 
interply layer. The maximum ayy observed in the +9 ply is less than 
10 percent of the applied stress. The a X y stress in the interply layer is 
negligible and can not be seen distinctly. Both a X y and c xx in the +9 
ply are seen to approach a finite value near the edge. A sign reversal is 
observed for these stresses. 

Interlaminar stresses q Y 7 , <77 7 . and a 7 y. - The interlaminar shear stress 
a xz has the greatest magnitude among these stresses. A steep increase in mag- 
nitude is seen near the free-edge zone. The stress magnitude in the interply 
layer is higher than that of the +9 ply. The interlaminar normal stress a zz 
appears to be of minor significance. However, it does have a definite trend. 
Two sign reversals are clearly noted. It starts as a compressive stress in the 
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interior of the ply or the interply layer, becomes a tensile stress for a short 
distance closer to the edge, and then reverts to a compressive stress near the 
edge. The maximum a zz observed is about 3 percent of the applied stress. 

Effect of width to thickness ratio on interlaminar stress peaks . - A lim- 
ited study with (±45) symmetric laminate has been conducted to assess the 
effect of width to thickness ratio on the interlaminar stress peaks. The 
stresses o xz and o zz are selected for the study because they appeared to be 
the only significant interlaminar stresses. The results are shown in figure 6. 
The stresses are computed 3/512 in. from the free edge. Four ratios of W/h, 
where W is the and h is the thickness, are chosen in the study. The 
results show significant sensitivity to the W/h ratio. Both stresses appear 
to have substantial magnitudes for the laminate with W/h = 8. As the lami- 
nate becomes thinner, however, o zz appears to reach a zero value while o xz 
continues to have a significant magnitude. 


Fracture Toughness Predictions of Composites for Individual 
and Combined Fracture of Modes I, II, III. 

Geometry and finite element idealization . - The three point bend test spe- 
cimens for end notch flexure (ENF) and mixed mode flexure (MMF) are used in the 
finite element evaluation of the fracture toughness for composites. These 
tests have been the subject of discussion and evaluation of ASTM D30.02 and 
D30.04 subcommittee meetings and specialty symposia sponsored by their subcom- 
mittees (ref. 3). Figure 7 describes the geometry and figure 8 shows the 
details of the finite element model with the crack-tip superelement. 

A line load was applied at the specimen midspan. The magnitude of this 
load was 120 lb. This magnitude corresponds to the experimental load which 
induced unstable crack propagation (ref. 3). The crack extension propagation 
was simulated by progressively deleting interply layer elements and repeating 
the FEA as described in reference 4. 

The composite material simulated was AS-graphite fiber/epoxy matrix (AS/E) 
unidirectional composite. The unidirectional ply properties, the interply 
layer thickness and its properties, and the appropriate material properties 
required for the three-dimensional finite element analysis were generated with 
the aid of ICAN (ref. 7). The properties for the Interply layers were assumed 
to be the same as those for the matrix. The specific values for the three- 
dimensional composite properties (using MSC/NASTRAN designation) and interply 
layer thickness predicted by ICAN are described in reference 3. 


RESULTS AND DISCUSSION 

The modes II SERR (strain energy release rate), using the local crack clo- 
sure method (ref. 3), is shown in figure 9 (dashed line) for 0.6 FVR. The 
curve determined using the global method is also plotted for comparison. The 
range of limited experimental data (ref. 3) is shown by dashed horizontal 
lines. Three points are worth noting in figure 9: 

(1) The global method predicts higher Gj j values than the local crack 
closure method. 
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(2) The two methods predict similar composite interlaminar fracture 
behavior. 

(3) A conservative assessment of the composite interlaminar fracture char- 
acteristics can be obtained by using the Gu determined from the local crack 
closure method. 

The mixed mode (I and II) composite interlaminar fracture characteristics, 
determined from the local crack closure method, are shown in figure 10 for a 
composite with 0.6 FVR. The fracture characteristics for each mode (mode I and 
11) are compared with those predicted for the mixed mode by using the global 
method (solid line) and the algebraic sum of the two modes from the local 
method (short dashed line). The following observations are worth noting in 
figure 10: 

(1) Both the local crack closure and the global methods predict similar 
mixed mode (mode I and II) fracture characteristics. 

(2) The local crack closure method predicts lower mixed mode G values 
than the global method. 

(3) Mode II fracture dominates the crack extension during slow and stable 
crack growth (a < 1.05 in.). 

(4) Mode I dominates the crack extension during unstable crack growth and 
rapid crack propagation (a > 1.15 in.). 

The general conclusions to be drawn from the above observations are: 

(1) The contribution of each fracture mode can be determined using the 
local crack closure method. 

(2) Mode I drives the composite interlaminar delamination in MMF 
specimens . 

(3) The local method provides a conservative assessment of mixed mode com- 
posite interlaminar fracture characteristics as was the case for mode II. 


Mixed Modes I, II, and III 
in Composite Laminates 

The computational procedure described in reference 3 has been modified to 
determine mixed modes I, II, and III in composite laminates (ref. 4). The lam- 
inate configurations used in these studies were unbalanced, and unsymmetric 
[-0 m /+0 n ]. These laminates were selected in order to evaluate the effects of 
the different laminate material-coupling coefficients, in the laminate force 
deformation relationships, on the individual and mixed mode fracture strain 
energy release rates (SERR). 

The specific laminates (-©36/+Q12 1 * 0 = 0°, 15°, 30°, 45°, 60°, 75°, and 
90°) were investigated for the laminate configuration effects on the SERR. At 
the 6 = 45° position, a total of three cases ( [ -45 m /45 1 2 ^ , m = 36, 60, and 84) 
were investigated for the interlaminar crack location effects on SERR. It is 
important to note that these specific laminate configurations and the 
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interlaminar crack locations were selected only for computational simulation 
convenience. They represent just one application of the present procedure. 


Effects of Ply Orientation 

The effects of ply orientation on the maximum individual and mixed mode 
fracture SERR (assuming a 480 lb load) are plotted in figure 11 for the seven 
cases of the C-©36 /+ ©123 <0 = 0°, 15\ 30°, 45°, 60° , 75°, and 90°) AS/E lam- 
inate. The relative dominance of the opening fracture mode SERR (Gi) and the 
negligible contribution of the tearing fracture mode SERR (Gm), on the maxi- 
mum mixed fracture mode SERR (Gy) are clearly observed in this figure. 

Another observation is that the "maximum" magnitude of mixed mode fracture (Gy) 
levels off at ply angle orientations greater than 60°. 

The above observations lead to the following conclusions: (1) the rapid 

or unstable interlaminar crack growth is dominated by the opening fracture 
mode; (2) the tearing fracture mode SERR (Gm) is negligible compared to Gj 
for ply orientation angles greater than 60°; (3) the individual fracture modes 
(Gi and Gn> and the mixed mode fracture are practically independent of ply 
orientation angle greater than 60°; and (4) ply orientation angles less than 
60° have significant influence on the SERR of individual fracture mode and 
mixed mode fracture. 


Effects of Laminate Configuration 

The effects of the laminate configuration on the SERR due to a 480 lb load 
are also shown in figure 11. The effects on the individual fracture mode 
SERR's are plotted in figure 11(a) to (c). The SERR for the mixed fracture 
mode (Gy) is plotted in figure 11(d). The important observations in this fig- 
ure are: (1) the shearing mode fracture SERR (Gn) appears to reach a maximum 

and then decrease with increasing crack opening for practically all the lami- 
nate configurations, (2) the tearing mode fracture SERR (Gm), on the other 
hand, continues to increase with crack length for some laminate configurations. 

The conclusion from the above discussion is that laminate configurations 
can be selected for "stable" shearing and tearing fracture mode crack growth 
for given composite components and loadings. It is important to keep in mind 
that advantages of this can be taken only in the absence of opening mode 
fracture. 


SUMMARY OF RESULTS 

The significant general and specific conclusions of the present investiga- 
tion are listed below: 


General Conclusions 

1. Finite element substructuring methods are computationally effective in 
dealing with the issues related to singular/near singular stress fields encoun- 
tered in advanced composite structures. 
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2. These methods are used at NASA Lewis Research Center to analyze a vari- 
ety of problems of current interest to the composites community. 

3. The results indicate the versatility of the finite element substructur- 
ing methods in the computational simulation of the complex composite singular- 
ity problems. 


Specific Conclusions 

1. The wave propagation velocities can be estimated from the early time 
displacement propagation response,, the estimates are in fair agreement with 
the theoretical predictions. 

2. The dynamic and static peak load stress contours are almost identical. 
The static peak load stress magnitudes are slightly higher. Stress predictions 
based upon a quasi-static approach lead to conservative estimates. 

3. The notch-tip region develops severe stress concentrations and any of 
the three stresses could cause or initiate a failure. 

4. The interlaminar stresses a xz , o zz , and o Z y are of second order 
compared to a xx in general. However, the interlaminar shear stress o xz 
could become substantial in the interply layer free-edge region depending on 
the configuration of the angleply laminate. 

5. Progressive mesh refinement is an effective way to describe the behav- 
ior of the interlaminar stresses as the free edge is approached. The results 
indicate stable convergence of these stresses with some noise as the free edge 
is approached. 

6. The magnitudes of the interlaminar stress peaks are sensitive to the 
width to thickness ratio of the laminate. Results obtained from relatively 
narrow specimens are not indicative of practical laminate behavior. 

7. The individual and mixed mode fracture SERR can be readily determined 
using a computational simulation procedure that consists of three-dimensional 
finite element analysis and integrated composite mechanics. 

8. Individual and mixed mode fracture SERR magnitude of [-0 m /+0 n 3 are 
strongly influenced by crack length, ply angle, and interlaminar crack loca- 
tion. However, the maximum magnitude of the mixed mode fracture SERR (Gx> is 
practically independent for ply orientations greater than 60° 

9. The tearing fracture mode SERR (Guj) has the smallest magnitude com- 
pared to opening (Gj) and shearing (Gn> fracture modes for this case. The 
tearing fracture mode is generally present in combinations with other fracture 
modes . 

10. Stress magnitudes ahead of the crack tip can be compared with corre- 
sponding local laminate strengths in order to determine the dominant stress 
which drives the crack. 


11. Collectively the results demonstrate that the procedure described 
herein can be used to computationally simulate/evaluate mixed mode fracture 
toughness parameters in composite components subjected to complex loadings. 
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FIGURE 1. - ASTM CHARPY TEST SPECIMEN AND FINITE ELEMENT IDEALIZATION 
OF SMOOTH AND NOTCHED SPECIMENS. 


DISPLACEMENT PROPAGATION 

3 MICROSECOND FRAMES 5 MICROSECOND FRAMES 



S/GLASS/EPOXY 

FIGURE 2. - DYNAMIC DISPLACEMENT PROPAGATION IN THE THREE COMPOSITE SYSTEMS SELECTED AFTER 3 ps AND 5 MS. 
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°xx °yy °*y 

FIGURE 3. - DYNAMIC/STATIC STRESS CONTOURS NEAR THE NOTCH REGION AT PEAK LOAD FOR S-GLASS/EPOXY NOTCHED FLEXURAL 
SPECIMEN. 
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- FREE-EDGE SUPERELEMENT 



FIGURE 4. - ANGLEPLIED LAMINATED PLATE GEOMETRY AND THE FINITE ELEMENT 
IDEALIZATION. 
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FIGURE 5. - FREE-EDGE STRESSES IN +10 u AND INTERPLY LAYERS OF (±10) LAMINATE. (SOLID CURVE, INTERPLY LAYER; DASHED 
CURVE, CNETER OF THE PLY). 




MIXED-MODE-FLEXURE (MMF) -MIXED MODE (I II) 

FIGURE 7. - SCHEMATIC OF FLEXURAL TEST FOR INTERLAMINAR FRACTURE MODE TOUGHNESS. 
NOTE: ORIGIN AT LEFT SUPPORT BOTTOM. 




FIGURE 8. - FINITE ELEMENT IDEALIZATION OF THE FLEXURE TEST SPECIMENS AND THE BACK-LIP SUPERELEMENT. 
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CRACK LENGTH, a, IN. 

FIGURE 9. - END-NOTCH FLEXURE ENERGY RELEASE RATE COMPARISON. 



GLOBAL MIXED (I & II) 
LOCAL MIXED (I & II) 
LOCAL I 



CRACK LENGTH, a, IN. 

FIGURE 10. - MIXED-MODE-FLEXURE ENERGY RELEASE RATE AND 
COMPONENTS (AS/E). LOCAL CLOSURE METHOD. 



FIGURE 11. - EFFECT OF LAMINATE CONFIGURATION ON STRAIN ENERGY RELEASE RATES ([-0 36 /+e 12 }; 
AS/E WITH 0.6 FVR). 
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simultaneously with the global response. The results from the select examples 
indicate that finite element substructuring methods are computationally 
effective for composite singularity mechanics. 
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